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1 Introduction 


This report discusses the application of a multigrid algorithm to the solution of the following one 
dimensional singular perturbation problem: 

- eu " + b(x)u l = /, 0 < x < 1 (1) 

with 

1 » e > 0, u(0) = uq, u( 1) = «i. 

Many other authors have discussed the application of various methods of solution to the alge- 
braic problem; in particular see Dorr [Dorr70a], Babuska [Babu69a], Ervin and Layton [Ervi85a], 
and Kellogg and Tsan [Kell78a]. 

Much of the literature regarding multigrid methods restricts itself to the solution of nice prob- 
lems. Indeed most authors require that the linear system be well conditioned in addition to symmet- 
ric and positive definite. However, in the case of singular perturbation problems as the coefficient of 
the second order term tends to zero the usual symmetric discretization fails to be of positive type. 
Thus the first measure taken in the numerical discretization of these problems is to replace the usual 
symmetric difference of the first order term with some form of skewed differencing. In particular 
for problems with turning points this may lead to a system of equations which is ill-conditioned 
for small e. 

From the standpoint of calculating a numerical approximation to the solution of problem (1) 
the first question is : does the discretization converge to the continuous solution? Then, assuming 
it does, how does the multigrid algorithm perform as a solver for this system of linear equations? 
What modifications, if any, are necessary in the multigrid algorithm? 

The main result shows that if the original system of equations resulting from the discretization 
of problem (1) is of positive type then the theoretical results for the multigrid algorithm developed 
in Kamowitz and Parter [Kamo85a] and in McCormick and Ruge [McCo82a] apply. 

From a computational perspective it is convenient to use for the coarse grid operators of the 
multigrid algorithm the operators that are the finite difference analogue of the original operator. 
In Section 6 the rate of convergence of the algorithm using the finite difference version of the coarse 
grid operators is considered. It is shown that the new rate of convergence is an 0(h 2 ) perturbation 
of the rate obtained using the Galerkin choice. 
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2 The Discrete Problem 


Three model problems are chosen for study, one where the sign of b(x) does not change and two 
where 6(x) changes sign. The three problems are designated 

Problem BL 

Lu = — eu 11 + u 1 = 0, 0 < x < 1 (2) 

Problem TP-1 

Lu = — eu n + {x — —)u ! = 0, 0 < x < 1 (3) 

Problem TP-2 

Lu = —eu" — (x — bu* = 0 0 < x < 1. (4) 

For all three problems the boundary conditions are 

u( 0) = 1, u(l) = 3. 

For the discrete problem as usual let N > 0 be chosen and set h = 1 /{N + 1). The interval 

n = (o, i) 


is discretized to form 

Q h = {ih : 1 < * < N}. 

The notation x,- refers to the point ih and u,' refers to u(x<). In addition the usual notation for 
finite differences is used: 


D+Ui = 


U,+ i - Ui 


D-Ui = 


Ui - Ui - 1 


DoUi = 


«,•+ 1 ~ u i - 1 


D+D-Ui = 


u,+i - 2 Ui + u,_i 


2 h ’ h* 

For problem (BL) the following two discretizations are considered: 


L\ui = -eD+D-Ui + D-Ui (upwind differencing) 


L\ui = 


—e 

T+hJ 2e 


D+D-Ui + D-Ui 


(see Kellogg and Tsan). 
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Since problems (TP-1) and (TP-2) have turning points at x = | it is necessary to change the 
direction of the discretization of the first order term. For problem (TP-1) the discretization tested 
is 


N 4- 1 

L z h Ui = -eD+D-Ui + (x, - \)D + u,, 1 < » < — - — (5) 

Ifa = - eD + D.u f + (x, - + 1 < « < N. (6) 

Similarly, for problem (TP- 2) the discretization is 

N + 1 

Liu, = -eD+D-Ui - (x, - \)D-u „ 1 < * < — - — (7) 

N 4- 1 

L\ui = -eD+D-Ui - (x,- - \)D + U{, — \-l <i < N. (8) 


Note that each of the discretizations Lj, k = results in a tridiagonal linear system of 

equations whose coefficient matrix may be denoted by 

L h = [ -Off Pi -lx ], 


with 

<*i > 0, 7 ,* > 0, and ft > + 7 1 *. 

Thus each A: = 1 , ... ,4 is an M-matrix and the linear system of equations 

Liu, = 1 <i<N, 


uq and ujsr+i fixed has a solution; see for example Berman and Plemmons [Berm79a]. 


3 Review of the Multigrid Algorithm 

The particular multigrid algorithm used to solve (1) has been discussed in detail in Kamowitz and 
Parter [Kamo85a]; thus only a cursory review is given here. In particular the details of the theory 
behind the convergence results can be found there. What is important to realize is that although 
the algorithm and convergence results discussed in the previous work apply to well conditioned two 
point boundary value problems the same bounds on the rate of convergence apply to the singular 
perturbation problems discussed here. 

In order to completely describe the algorithm a number of spaces and operators need to be 
defined. For the various spaces choose g nested grids 

n P h 
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where 


p = 2 9 " 1 . 

Denote by Skh the space of grid functions defined on Qkh- From now on the notation kh will be 
replaced by k. Denote by Gk the smoother, the restriction map, and the interpolation 
map. Associated with each space S k there is also a nonsingular tridiagonal operator 


Lk : S k S k> 


where L k is of positive type. The operators L k are again denoted 

L k = [ -ft ft -ft ). 

These operators will be formally defined after the statement of the algorithm. 

The following is an outline of the algorithm used. Assuming U is the nth approximation to 
the solution of the system of equations 

LiU x = F X 


algorithm MG(U£, L k) F k ,k) returns the next iteration of the multigrid algorithm. The grid 

layer is denoted by fc; set A: = 1 to start the algorithm on the finest grid layer. 

Algorithm MG(U£, L ky F k , k) 

1. Coarse Solve: If k = g (coarsest grid) then return 

MG <- Lg-'Fg 

otherwise 

2. Smooth: Apply the smoother, G *, call the result of this step U£. 

If k = 1 use Ui as the initial guess, otherwise use 0. 

3. Recursively Apply the Algorithm: Set 

u n+1 <- u? + I k k+1 MG( 0, L k+ u I k k +1 (Fk - L k U£),k + 1). 

4. (optional) Smooth Again: Set U n+1 <- G k U n+1 . 

5. Return: Set MG *- U n+1 . 

As defined algorithm MG is called a symmetric V-cycle if step 4 is used; otherwise algorithm 
MG is referred to as a slash cycle (following the notation of McCormick and Ruge) . 

For the smoother, G, choose m-applications of damped Jacobi iteration with parameter a. 
Formally, repeat for 1 < r < m 


tr - v;+ - w% 

The interpolation operator, J* +1 , is defined as follows. For points common to Qk+i and to f2* 
set 

!4‘ + it'W = v, 
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and at odd (new) points of fi*. require 

{Lk[lUiU)} 2j _ x = 0. 

This results in the explicit system of equations at the point 


Tjk 

U 2j-1 



The restriction operator, /£ +1 , is 



.^23 -1 


Tjk 

U 2j- 


1 + U 2j + 


- 2l -Cf, fc 


0k 2,+l 

^2 3+1 J 


Note that if L t is symmetric then 



A fundamental observation due to McCormick and Ruge [McCo82a] in the symmetric case is that 
S'* can be written 

Sfc = Range I* +1 © Nullspace (9) 

For n on-symmetric problems the above decomposition follows directly from the characterization of 
Range J* +1 and Nullspace t; see Kamowitz and Parter. 

Finally the coarse grid operators, Lk+i, are chosen by setting 


L k+1 = L k+1 = I k k +1 L k I k k+x . 


A direct computation shows that 


£jk+l — 


- a 


,fc+i 


3 


-r 


,Jfe+l 


with 


a k+1 

i 


ok+1 

3 


-rr 


44'- 1 


4-1 


4 - 


44-i 


4-i 


44+1 

4 +i J 


44+i 


4 


2;'+l 


( 10 ) 

( 11 ) 

( 12 ) 


Note that the choice of is the ‘natural’ choice; however it is not the only feasible choice for 
£*+i- 
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3,1 Convergence Theory 

For completeness some results from [Kamo85a] are repeated. 

Let 

e n - U true - U n 

be the error in the nth iterate. Here U trw , corresponds to the true solution of the algebraic problem 

LfrUtrue — /• 

The rate of improvement of the nth iterate of algorithm MG is then 


cn-II 


To estimate p , the asymptotic rate as n — ► oo, one needs to bound the p n . 

For later use define 

11 * 111 * = (L h x,x) = Y {Lhx)^. 

i 

First the two grid process is considered. Given an initial error 

6 ° = u true - U°, 

the two grid process yields: 

1. Smoothing 

0 ~0 0 
6 — > 6 = G € 

where G* corresponds to the linear part of the smoother G. Note that from (9) 

~e° = r] + I% h w 2h 

where 

i ) G Nullspace I% h Lh and W 2 h € 

2. Restricting the residual r* = L*e° to yields 

R2h = ll k L h e° = ll h L h r) + I 2 h h L h l£ h w 2h = I 2 h h L h l£ h w 2h 


since 


T) G Nullspace 1*^ L* 
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3. Computing the coarse grid correction directly (as is the case for the two grid algorithm) is 
equivalent to solving 

L 2 h'/>2h = Rih — LihWih- 

Here 

L 2h = Uh = iTlA 

SO 

1p2h = U>2 h- 

4. Finally, correcting ?° using the coarse grid correction yields e 1 : 

e 1 = Ut^-U 1 

= U true -(U 0 + I^2h) 

= e° - ihfoh 
= n + I 2 h w 2h ~ l2h^2h 
= V- 


Note that if 

6° G S 2h 


then one step of the two grid procedure solves the problem! In general, since ||G'|| < 1 the rate p l 
satisfies 




In this case 



so 



Notice that the r) term is related to the action of the smoother while the I 2 h w 2h term is eliminated 
by the multigrid process itself. 

In Kamowitz and Parter [Kamo85a] an explicit decomposition of Sh was found in terms of the 
eigenvalues and eigenvectors of the damped Jacobi scheme. This decomposition was exploited to 
compute bounds on the rate of convergence. 

For the two grid scheme, the bounds on p> for a given m and a, are given in Table 1. Note that 
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I 


m 

a 

.333 

.500 

.667 

.750 

1.00 

1.333 

i 

.500 

.333 

.400 

.429 

.500 

.572 

2 

.250 

.111 

.160 

.184 

.250 

.326 

3 

.125 

.078 

.088 

.093 

.125 

.187 

4 

.068 

.062 

.068 

.072 

.083 

.109 


Table 1: Predicted Rates for 2 Grids 


m 

a 


.333 

.500 

.667 

.750 

1.00 

1.333 

i 

.633 

.577 

.561 

.561 

.577 

.614 

2 

.435 

.408 

.417 

.424 

.447 

.475 

3 

.336 

.335 

.349 

.357 

.378 

.403 

4 

.283 

.293 

.307 

.314 

.333 

.357 


Table 2: Predicted Asymptotic Rate 


.5 case will 
[McCo82a] 

4 The Singular Perturbation Case 

For the study of singular perturbation problems it is necessary that if L ^ is an M-matrix then the 
coarse grid operators & = 1,2, .. . are also M-matrices. This is necessary to insure that each 

subproblem 

£*+ 1^+1 = fk+i 

is solvable. In particular since the smoother Gfc+j depends on £jb+i, then if Ljb+i is an M-matrix, 
!l<3fc+i|| < 1. 

Lemma 4.1 If 

Lk = [ -oj % 1j ) 

is an M-Matrix then 

L k+1 = [ — a* +1 /?J +1 -7? +1 ] 

is also an M-matrix, where a* +1 , /3j +1 and 7y +1 are given by equations (10-12). 


the optimal rate is obtained for a = .5. In the succeeding sections the bound in the a = 
be obtained experimentally for problems ( BL) y (TP-1), and (TP-2). 

In the case where the number of grids is arbitrary, estimates based on the ideas in 
are given in Table 2. As indicated in [Kamo85a] these bounds are not sharp. 
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Proof: The hypotheses that Lk is an M-matrix insures that 


a’j > 0 , > 0 , and 7 * > 0 

so aj +1 and 7 y +1 are also positive. In addition, /?y +1 must satisfy 

(3j +1 > aj +1 + 7 y +1 . 

Recalling the expressions for and 7 y note that this is equivalent to requiring 



44-1 

4-1 


44+1 > 44'- 1 44+ 1 

4 + i " 4 + i 


or, 



+ 


4 

P 2j+1 



which follows directly from the fact that L* is diagonally dominant. 

Since each of the Ljb+i is an M-matrix, the theory developed in Section 3.1 can be applied in 
the singular perturbation case. 


5 Experimental Results 

The algorithm of the previous sections was applied to problems ( BL) y (TP-1) and (TP-5). Prob- 
lem \BL) exhibits a boundary layer at x — 1. From a computational point of view the system of 
linear equations that is solved is well conditioned even for small e\ however the fact that the linear 
system is not symmetric leads to computational difficulties in computing the experimental rate of 
convergence. 

In the case of problem ( TP-1) there are two boundary layers; at x — 0 and at x = 1. In addition 
the system of equations that is solved is ill-conditioned. 

The discrete problem corresponding to problem (TP- 2) is well conditioned. The fact that there 
is an interior turning point at x ~ 1/2 does not appear to lead to computational difficulties. 

By using the one sided difference schemes discussed in the previous section the linear systems 
arising from the discretization of the problems are of positive type. 

From a practical standpoint, however, what effect does the ill-conditioning have on the observed 
solution? In particular, how should the rate of convergence (reduction in the error) be measured? 
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5.1 General Remarks Regarding the Experiments 

In all the experiments N, the number of points on the fine grid, was 63. The initial guess, Z7°, was 
constructed so that the initial error could be chosen in advance. In other words, 

Uj = «(*;) - e r 


For the experiments discussed here 

ey = (-I)'. 

The number of smoothing steps for the experiments discussed here was set to one. Experi- 
ments were run where the number of smoothing steps was greater than one. When more than one 
smoothing step was used there was no observed qualitative difference in the behavior of the algo- 
rithm compared to the behavior for non-singular perturbation problems. Unless otherwise noted 
the results are for the two grid case. The computer used was a VAX-11/780 and the tests were run 
in double precision (roughly 16 decimal digits of accuracy). 

The damped Jacobi parameter, a, was chosen to be .5. The following heuristic argument due 
to Brandt [Bran77a] suggests why a — .5 is optimal for one smoothing step. His suggestion is 
to choose a so that the range of eigenvalues which are reduced by the damped Jacobi process is 
equal to the range which is left alone by the process. As noted in [Kamo85a] the eigenvalues of the 
damped Jacobi scheme come in pairs A(/i) and A(/i) where 


A(m) = 


a 

1 + a 


and 


A(jx) = 


a — fi 
1 + cl 


and fi are the eigenvalues of the scheme for a = 0. The eigenvalues fi are real and distinct and as 
— > 0 they fill out the interval [0,1]. Brandt’s requirement corresponds to choosing a so that 


-A(-l) = A(0) 


or in other words to taking a — 1/2. Note that this is equivalent to requiring 

-A(l) = A(0). 

Indeed in [Kamo85a] for 1 smoothing step the theoretical and experimental results indicate that 
this choice of a is optimal. 


5.2 Results for the Boundary Layer Problem 

For problem ( BL ) the solution to the analytic problem can be calculated explicitly. For all e > 0 
the solution u £ (x) is 

u e (x) = C\e x l c + C%. 
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E 

Grid 

l 

.1 

.01 

.005 

.003 

.001 

0 


.62 xlO - ^ 

.11 xlO -2 

.54 xlO -2 

.73 xlO' 2 ^ 

.86 XlO' 2 

.011 

.012 


.25 Xl0 -2 

.45 xlO -2 

.18 XlO' 1 

.21 xlO' 1 

.28 xlO' 1 

.027 

.024 


.98 Xl0 -2 

.12 xlO -1 

.37 XlO' 1 

.44 xlO' 1 

.46 XlO' 1 

.048 

.048 


.39 xl0 _1 

.41 xlO' 1 

.80 xlO' 1 

.92 xlO' 1 

.96 xlO' 1 


.1 

5 

.15 

.15 

.19 

.20 

.21 

.21 

.21 

6 

1 

1 

1 

1 

1 

1 

1 


Table 3: Condition Number - Problem BL 


The constants Ci and C 2 satisfy 


C x = 



C 2 = 1 - Cl 


e 1 /' - 3 
e 1 /* - 1 ‘ 


Note that 


lim Ci = 0, lim C*i — 1 

0 e-0 


SO 


limti c (x) ~ 1. 
e — ►O 


Denote by Uh )€ (x) the numerical solution of problem (BL) for a fixed h and e. The consistency 
condition on the solution then requires that for fixed e 


lim U hte (x) = u e (x). 
/i— ►o 


Dorr [Dorr70a] proved that for fixed h 


YimU h>e (x) = 1. 

The discretization of problem ( BL ) using L\ (standard upwind differencing) results in an ap- 
proximation to the true solution which is an 0(h) approximation. To improve on this approximation 
Kellogg and Tsan point out that the use of L\ results in an 0(h?) approximation for e > 0. For 
the reduced problem (e = 0) L\ gives an 0(h) approximation. 

It is important to note that the linear system of equations arising from the discretization of 
Problem [BL) is not ill-conditioned. Table 3 displays the LINPACK [Dong79a] estimate RCOND 
(an estimate of the inverse of the condition number) for L \ , h = 1/64 and e = 1, .01, .005, .003, 
and .001. 

First, some general conclusions about the experimental results. For both L\ and L\ the al- 
gorithm converged to the solution of the algebraic problem with an asymptotic convergence rate 
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identical to the rate predicted in Section 3.1. Moreover there was no observed qualitative difference 
in the behavior of the algorithm with respect to the choice of L £ or L \ . However, as e — > 0 the 
behavior of the iterates changes dramatically. 

Define 


i?n — 


INI 

lkn-1 


II 


where 


r n = F - L h U n 


is the residual after the nth iteration and the norm used is the I 2 norm. From a computational 
point of view this is a convenient measure of the rate of improvement since one does not know the 
true solution. From past computational experience this ratio is bounded above by the theoretical 
error reduction rate. Indeed for large e, say e — 1, this is the case. For small e, say e = .001, after 
a small number of iterations R n exceeded the rate predicted in Table 1, and then as n increased R n 
declined towards the predicted rate; see the solid line in Figure 1. The dashed line will be referred 
to later. In both cases two grids are used. 

An explanation for this behavior is that as the algorithm proceeds the fact that the error is 
being measured in an asymmetric norm causes problems. The user of the algorithm should be 
careful to note that while in principle all norms on finite dimensional spaces are equivalent the use 
of a symmetric norm results in much better observed behavior of this particular algorithm. 

To demonstrate this hypothesis the problem 


L h U = F 


is transformed into the equivalent symmetric problem 

D~ 1 L h D(D~ l U) = D _1 F = F. (13) 

The matrix D is a positive diagonal matrix whose entries are given by 


d\ = 1, 



Applying this transformation to L* results in L s h which is now symmetric. The matrix L’ h is 
denoted 

L s h =[ -aidi-i/di Pi -7id»+iM ]• 

Denote by algorithm MG S algorithm MG applied to problem (13). The error resulting from 
applying algorithm MG S is related to the error from applying algorithm MG as follows: 
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0.10 


0.00 


T— I 1 I 1 V 1 1 1 1 I I I I 1 1 1 1 I 

1 6 11 16 21 26 31 36 41 46 51 56 61 66 71 76 81 86 91 

Iteration 


Figure 1: J? n , for e = .001 
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Lemma 5.1 If e n is the error after the nth iteration of applying algorithm MG then the new error , 
e n+1 ; is related to the error in the symmetric problem by 

(MG s )(D~ 1 e) = D~ l MGe = D~ x e n+l . 


Proof: By the definition of e n , 


and 


e n = U -U n 


e n+1 = U - U n+1 = MGe n . 

Also, 

D~ x e n = D~ X U - D^ 1 U n , 

Applying each step of algorithm MG S to D~ l U n results in D~ 1 U n + x as a straightforward calcula- 
tion shows. Hence 

MG 3 D~ l U n = D~ l MG U n . 


Applying this transformation to the error e n and computing the rates 


D-'Rn = 


I \D-^\\ Lk 
Wd-^-'Wl 


results in the more usual behavior displayed by the dashed line in Figure 1. Figure 2 displays the 
location where the maximum norm of the error is taken on versus the iteration number. Notice 
how as the iteration proceeds the location where the maximum occurs ‘drifts. ’ Applying D~ x and 
then measuring R n has the effect of ‘fixing’ the location where the maximum error is taken on. 

Unfortunately the entries of the transformation matrix D satisfy 


Di — ► oo 

as e — ► 0 so this is not a stable method for solving this problem in general. 

In summary, for problem (2) where the coefficient matrix was not poorly conditioned but was 
non-symmetric the only computational difficulty was in choosing the norm in which to measure 
the rate of improvement of the algorithm. 


5.3 Turning point Problem - 1 

In this section consider the solution of 


— £u e n + (x — 



(14) 
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Figure 2: Index of Maximum Error, Problem BL 
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€ 

Grid 

l 

.1 

.01 

.005 

.003 

.001 

i 

.00023 

.00032 

.20 Xl0~ 6 

.19 xl0-^“ 

.13 xl0“ 12 

.71 xlO' 21 

2 

.0023 

.0013 

.76 Xl0~ 6 

.66 xlO" 9 

.41 xlO" 12 

.23 XlO' 18 

3 

.0096 

.0085 

.83 Xl0~ 2 

.85 xlO -2 

.80 XlO" 2 

.60 XlO' 2 

4 

.038 

.037 

.46 XlO -1 

.55 xlO -1 

.60 XlO" 1 

.65 xlO -1 

5 

.15 

.146 

.16 

.17 

.18 

.19 

6 

1 

1 

1 

1 

1 

1 


Table 4: Condition Number - Problem TP-1 


with 

u e (0) = 1, u £ (l) = 3. 

This problem has two boundary layers, at x — 0 and at x = 1. The asymptotic solution satisfies 
(see [Krei74a] and [Dorr70a] ) 

lim u £ (x) = 2. 

For the discrete problem again 

L s h U = F (15) 

is solved where L\ is given by equations (5) and (6). The condition number of L\ tends to infinity as 
e — ► 0 for a fixed h. Table 4 displays the LINPACK estimate RCOND for h— 1/64 and e = 1, 
.01, .005, .003, and .001. It is important to note that although the operator corresponding to the 
original discretization on the /i— grid is ill-conditioned, Ljt, k = 3, 4, . . . are not ill-conditioned. In 
fact in the limiting case where the coarse grid has one point the system of equations to be solved is 
a 1 X 1 system which has condition number 1. This feature of the multigrid algorithm is reassuring 
to the user since it means that the actual system being solved directly in step 1 of algorithm MG 
is well conditioned and no special measures need to be taken. 

Although the coarse grids used in the algorithm do not resolve the boundary layers the algorithm 
still converges. This is because the role of the coarse grids is to solve for the error in the solution, 
rather than to represent the solution itself. The boundary layers are not seen in the expression for 
the error. The observed rate of convergence was independent of the value of e used. 

For all tested values of e the multigrid algorithm converged to the solution computed by Gaus- 
sian Elimination to the algebraic problem. However, the ill-conditioning of the algebraic system 
resulting from the discretization of the continuous problem opens up the question of to what so- 
lution does the multigrid algorithm converge? Since the original system is ill-conditioned there is 
a family of solutions U for which the residual is small. Indeed for e = .001 the condition number 
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Grid 

i 

.1 

.01 

.005 

.003 

.001 

0 

1 

.64 XlO -3 

.98 XlO -3 

.45 xl0~ 2 

.69 xl0~ 2 

.89 xlO” 2 

.013 

0 

2 

.26 xlO -2 

.39 xlO -2 

.17 xlO" 1 

.24 xlO" 1 

.28 xlO -1 

.031 

0 

3 

.99 xicr 2 

.11 XlO- 1 

.25 xlO’ 1 

.33 xlO -1 

.37 xlO" 1 

.038 

.035 

4 


.040 

.063 

.080 

.089 

.095 

.096 

5 

.15 

.15 

.18 

.20 

.21 

.22 

.22 

6 

1 

1 

1 

1 

1 

1 

1 


Table 5: Condition Number - Problem TP-2 


estimate is .71 xlCT 21 which is less than the unit roundoff of the VAX-11/780. Thus there is no 
reason to expect that the solution calculated by solving equation (15) is an adequate approximation 
to the solution of the continuous problem. In the experimental runs that was indeed the case and 
both algorithm MG and Gaussian Elimination returned 0 for the solution. 

Since for this problem the error is not skewed as in problem ( BL ), it is not surprising to see 
that the ratio of improvements R n are bounded above by the predicted rate from Table 1. This 
can be seen in Figure 3, for e = .001. 

5.4 Turning point Problem - 2 

In this section the solution to 

- eu" - (x - l)u £ ' = 0, (16) 

is discussed. The asymptotic solution, see [Krei74aj is 

u e (x) - 1, 0 < x < i 

u ff (s) = 3, ^ < x < 1. 

Although one sided differencing, in this case is used to compute the solution, the system of 
equations that is solved is not ill-conditioned. As in Section 5.3 the LINPACK estimate RCOND 
of the condition number was computed; see Table 5. 

It appears that the condition number is related to the boundary data. For problem ( TP-1 ) the 
discretization uses information from the boundary layer to estimate ti c (x) at interior points. The 
boundary data for problem [TP-2) is well represented (no boundary layers) in the interior. 

From the standpoint of the multigrid algorithm again as the grids get coarser the condition 
number improves. Also, as for problem (TP-1) the rate of convergence of the algorithm was 
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Figure 3 : Rn~ Problem TP-1, e = .001 
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independent of the value of e, and was identical to the rate predicted in Section 3.1. In addition 
R n was bounded above by the predicted rate. 


6 Comparison of I^Lhlih an d L 2 h 


From a computational point of view one would like to choose for the coarse grid operator L^h the 
tridiagonal operator analogous to L ^ obtained by finite differences. This choice for L 2 /i will be 
denoted L For now only the symmetric problem 


- (pu'y = /, p(x) > po > o 

with boundary conditions 

u(0) = u(l) = 0 

is considered. Here p(x) is a smooth function. 

Lemma 6.1 The operator L ££ is related to L^h by 

L 2h = L f 2 d h + h 2 [-A j By —Cj ] 
where Aj, Cj depend on p(x), p'(x) and p"(x), and Bj = Aj + Cj. 

Proof: Recall that the coefficients of L^h are given in equations (10-12). In particular 


(17) 


— 


«2t«2fc-l 


02k- 


1 J 


Expanding 


. 4 & — i . 

<*2 Jfc = P( 7, h )> 


«2k-l = ?( 


2 

4fc — 3 


h) 


and 


„ . 4fc — 1 . ,4fc-3 n 

02k— l = P( — h)+p{— — h) 

in Taylor series about the point x = (4A; — l)/i and collecting terms yields 


“* = o 


1 Ot2k<X2k-l 


2 02k— 1 
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_ i [p(«) + rp'(g) + ~ fp'C*) + jyjx) - c 2 ft 3 ] 

2 [p(*) + fP'(x) + T-p"( £ ) + ^l^ 3 ] + b(*) - | P f (z) + *TP"{x) - Cl * 3 ] 

_ 1 p(x) 2 + /i 2 [|p"(x)p(x) - ^(p'(x)) 2 ] + Dh t 3 

2 2p(x) + £p"(i) + Eh 3 

= ^p(*) + 

= Oi^ d + h 2 A k 

where A* depends on p{x) y p\x) and p"(x). 

A similar calculation yields 

% = 7 j l d + h 2 C k . 

In addition 

Pk = Oik + Ik 

= + h 2 A k + l[ d + h 2 C k 

= p{ d + h 2 B k . 

The error analysis in Section 3.1 shows that if L 2 h is not used then the coarse grid correction 
t/>2 h w iU not equal w 2 f l . Thus a bound on the eigenvalues of the eigenvalue problem 

A lit* = W 

is needed. 

Lemma 6.2 The eigenvalues A of 

AL& = W 

satisfy 

min 1 7, • | 

Proof: By Lemma 6.1 

L 2h = L{i + h 2 [ -Ay By -C-y], 

thus 

(A - 1 = /i 2 ([ -A,- B, -C7.- ]^>. 


20 



Summation by parts yields 


A — 1 = h 


= h 


< h 


2 ([ -Aj Bj —Cj ]M) 
2 ’LCijrpi - & + i) 2 
2 maxjCil 




min 1 7 / 

= Kh 2 

< OO 

because of the smoothness of p and since \^{ d \ > po > 0. In addition, since L^ d is positive definite, 

(A - 1 )(L f 2 ^,rp) > 0 


which implies 


Therefore 


Recall from equation (9) 
and if L 2 h is used then 
In the case where is used to solve 

the error is denoted e 1 . Now 


A - 1 > -Kh 2 . 
|A - 1 | < Kh 2 . 
e° = 7 + lh w 2h 


* = n. 


/ J A 

L 2h ^2h = L 2h w 2h 


e 1 = r) + I% h (w 2h - xp 2h ). 

Since Range I% h and Nullspace I^Lh are Lh — orthogonal (see [Kamo85a] ), 


l, ? 1 iiL = n»?iiL+ii^(^-^)ii2 k . 


21 


Applying Lemmas 6.1 and 6.2, 


ll^hO^h - V>2fc)llL = (Lhl2h(. w 2h - ^ 2 h), l2h( w 2h - V»2fc)) 

= 7£{L2h(u>2h ~ 'P2h),(*»2h. ~ ^2h)) 

= ^(^2fc(^2 h-(Li£) L2h^2h),( w 2h ~ (L{*) L2 h W2h)) 

= \{FL X £rv2 h ,FL^W2 h ) 


where 


F = (I-Ll{ 2 (L'i)- L\' 2 ). 

The eigenvalues of L\tf (L{f) l\£ are the same as the eigenvalues of (L^h) £ 2/1 > so by Lemma 6.2 




< I (Kh?)\L'£wi K ,L'£wn) 
= i(^fe 2 ) ! <£2hUl2h,W ! fc) 


Therefore 


and 


||&(“» - Mill, < («-A 2 ) 2 i|/ 2 \<»2»ii^ 


Wli, < llilli. + (*-a 2 ) 2 ii^u I ,»ii 2 Ii . 


Since ||G|| ifc < 1, 

Finally, the rate using is 


J^lk<ll*k + **Vk- 


fd _ H £ II Lk 

" - Mk 


< 


Lh 

k 


+ Kh 2 


Recall the rate using L, 2 h is 


so 

p fd <p+Kh 2 . 

Remark: Although the discussion in the previous section is restricted to symmetric problems, the 


same result extends to the general two point boundary value problem 


1 - (pt/y + bu 1 + 

cu = /, p(x) > po > 0, c(x) > 0 

(18) 

since problem (18) can be rewritten 

- (pti')' + cti = ? 

(19) 

where 

p = 

(20) 

1 

• ■ •<! 

(21) 

i 

| 

' ■ •<•’» 

(22) 

! 

The function p(x) is 



\ ?(*) = 



j Note that p, c and f are defined since 

p{x) > po > 0. 



6.1 An Illustration 

A > J 

As an illustration of the effect of replacing L 2 h with L J 2h consider using for the smoother G in 
algorithm MG one sweep of odd Gauss-Seidel smoothing. In other words 

{GU)i = Ui for t = 0 mod 2 

and for * = 1 mod 2 

{GU\ = + ~tiU i+ 1 + /,*]. 

Pi 

This guarantees that the error e lies completely in Rangel ^ In other words, the rj term is 0 since 
for i = 0 mod 2 

(Lkl2h u, 2h)i 0 
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yet 


0 = (L h T),I^ h W2h) 

= {ruLhlih^h)- 

In addition for points x t , with t = 1 mod 2, 

0 = r, = Lh* — (Lkv)i + (Lhl 2 h w 2 h)i 

= (LkVh 

which implies = 0 since r? t _x = *7i + 1 = 0. 

Thus 

£ = l2 h W 2 h 

for some If £ 2/1 is used then one iteration of the algorithm results in completely solving the 

problem. 

For a numerical example problem (17) is solved with 

p(x) = «*“(**). 

The right hand side / is chosen so that the the true solution is 

turtle (*) = sin(jrz). 

As expected when L^h is used the algorithm converges in one step. When is used then 
the observed rate of convergence corresponds to the Kh? term of Section 6. Table 6 displays the 
observed rate of convergence versus h for four different choices for the number of grids. The column 
p(h) corresponds to the observed rate of convergence for each value of h. The value of 5(h) is 


p{h) ‘ 

Since the error in the two grid case is 0(h 2 ) one expects that 8(h) — ► 4 as h — ► 0 for two grids. 
Indeed this is the case. The reason why p(h) varies with the number of grids used is that there are 
‘pollution’ effects from not solving each coarse grid equation exactly. In other words 

e 1 = Vh + I-ih r l2h + • • • + 0(h 2 ), 

where the 0(h 2 ) term corresponds to the error made by not using L^h- 
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2 grids 

3 grids 

4 grids 

5 grids 

h 

p0 0 

6(h) 

fW 

8(h) 

pW 

8(h) 

P(h) 

6(h) 

1/16 

.50 x 10' 2 


.27 x 10 _1 


.85 X 10" 1 




1/32 

.12 x 10" 2 

4.17 

.62 x 10“* 

4.35 

.28 X 10" 1 

3.04 

.86 x 10 _1 


1/64 

.30 x 10~ 3 

4.00 

.15 x 10"* 

4.13 

.65 X 10“ 2 

4.31 

.29 x 10 _1 

3.01 

1/127 

.75 x 10“ 4 

4.00 

.38 x 10" 3 

3.95 

.16 X 10“* 

4.06 

.66 x IQ” 2 

4.33 


Table 6: Rate using Lj* 
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